rm(list=ls())

pacman::p_load(lubridate, tidyverse, haven, gtrendsR, rdrobust, rddtools,
               broom, ggpubr, gridExtra)
select <- dplyr::select

setwd('put_your_wd_here')

# shootings
allshootings <- read_csv('datasets/list_of_shootings.csv',col_select = 1:4)
allshootings$date  <- mdy(allshootings$date)
top10 <- allshootings$shooting[allshootings$top10==1]

# Frontline Documentaries ----
df <- read_csv('datasets/frontline2.csv')
df$date <- ymd(df$date)
mindate <- min(df$date) + months(2)
maxdate <- max(df$date) - months(2)
dates <- allshootings %>%
  filter(date >= mindate & date <= maxdate & include==1)

# merge in guardian documentary
gua <- read_csv('datasets/gun_nation.csv') %>%
  dplyr::select(date, gun_nation)
df <- left_join(df, gua)

# GUN STORE ----
out <- readxl::read_excel('datasets/gun_store.xlsx', sheet=2)
out <- filter(out, Geography == 'US')
out$date <- ymd(out$Date)
out <- out %>% rename(gun_store = Views) %>%
  select(date, gun_store)
df <- full_join(df, out, by='date')

makeEstRDT <- function(out, treat, name, org, y, time=2){

  min <- treat - months(time)
  max <- treat + months(time) 
  out$running_day <- as.numeric(out$date - treat)
  trim <- out %>%
    filter(date < max & date > min)
  y = trim %>% select(!!enquo(y)) %>% pull
  
  house_rdd <- rdd_data(y=y, x=trim$running_day, cutpoint=0)
  bw_ik <- rdd_bw_ik(house_rdd)
  reg_nonpara <- rdd_reg_np(rdd_object=house_rdd, bw=bw_ik)
  fullest <- data.frame(Coef=reg_nonpara$coefficients,
                           lower=(reg_nonpara$coefficients - 1.96 * reg_nonpara$coefMat[2]),
                           upper=(reg_nonpara$coefficients + 1.96 * reg_nonpara$coefMat[2]),
                           shooting=name,
                           outcome=org)
  return(fullest)
}

# this function estimates effects in SDs for meta analysis
# makeEstRDT2 <- function(out, treat, name, org, y, time=2){
#   
#   min <- treat - months(time)
#   max <- treat + months(time) 
#   out$running_day <- as.numeric(out$date - treat)
#   trim <- out %>%
#     filter(date < max & date > min)
# 
#   y = trim %>% select(!!enquo(y)) %>% pull
#   
#   sd_est <- sd(y, na.rm=T)
#   house_rdd <- rdd_data(y=y, x=trim$running_day, cutpoint=0)
#   bw_ik <- rdd_bw_ik(house_rdd)
#   reg_nonpara <- rdd_reg_np(rdd_object=house_rdd, bw=bw_ik)
#   fullest <- data.frame(Coef=reg_nonpara$coefficients/sd_est,
#                         lower=(reg_nonpara$coefficients - 1.96 * reg_nonpara$coefMat[2])/sd_est,
#                         upper=(reg_nonpara$coefficients + 1.96 * reg_nonpara$coefMat[2])/sd_est,
#                         shooting=name,
#                         outcome=org)
#   return(fullest)
# }

gunned_down <- vector('list', nrow(dates))
for(i in 1:nrow(dates)){
  gunned_down[[i]] <- makeEstRDT(df, treat=dates$date[i],name = dates$shooting[i], org= 'Gunned Down', y=gunned_web)
}

nra_web <- vector('list', nrow(dates))
for(i in 1:nrow(dates)){
  nra_web[[i]] <- tryCatch({makeEstRDT(df, treat=dates$date[i],name = dates$shooting[i], org= 'NRA Under Fire', y=nra_web)},
                           error=function(e){NA})
}

gun_nation <- vector('list', nrow(dates))
for(i in 1:nrow(dates)){
  gun_nation[[i]] <- tryCatch({makeEstRDT(df, treat=dates$date[i],name = dates$shooting[i], org= 'Gun Nation', y=gun_nation)},
                           error=function(e){NA})
}

gun_store <- vector('list', nrow(dates))
for(i in 1:nrow(dates)){
  gun_store[[i]] <- tryCatch({makeEstRDT(df, treat=dates$date[i],name = dates$shooting[i], org= 'Gun Store', y=gun_store)},
                              error=function(e){NA})
}

dat <- bind_rows(
  bind_rows(gunned_down%>% discard(is.na(.))),
  bind_rows(nra_web %>% discard(is.na(.))), 
  bind_rows(gun_nation%>% discard(is.na(.))), 
  bind_rows(gun_store%>% discard(is.na(.)))
) %>% na.omit() %>%
  remove_rownames()

# estimate meta analytic estimates
# dat4meta <- bind_rows(
#   
#   makeEstRDT2(out, treat=dates[2],name = 'Orlando',org= 'Gunned Down', y=gunned_web),
#   makeEstRDT2(out, treat=ymd('2015-10-1'),name = 'Umpqua',org= 'Gunned Down', y=gunned_web),
#   makeEstRDT2(out, treat=ymd('2015-12-2'),name = 'San Bernardino',org= 'Gunned Down', y=gunned_web,1),
#   
#   makeEstRDT2(out, treat=dates[3],name = 'Stoneman Douglas',org= 'Gunned Down', y=gunned_web),
#   makeEstRDT2(out, treat=dates[4],name = 'Atlanta + Boulder',org= 'Gunned Down', y=gunned_web),
#   makeEstRDT2(out, treat=dates[5],name = 'Gilroy + El Paso + Dayton',org= 'Gunned Down', y=gunned_web),
#   makeEstRDT2(out, treat=dates[6],name = 'Las Vegas',org= 'Gunned Down', y=gunned_web),
#   
#   makeEstRDT2(out, treat=dates[4],name = 'Atlanta + Boulder',org= 'NRA Under Fire', y=nra_web),
#   
#   makeEstRDT2(out, treat=dates[3],name = 'Stoneman Douglas',org= 'Gun Nation', y=gun_nation),
#   makeEstRDT2(out, treat=dates[4],name = 'Atlanta + Boulder',org= 'Gun Nation', y=gun_nation),
#   makeEstRDT2(out, treat=dates[5],name = 'Gilroy + El Paso + Dayton',org= 'Gun Nation', y=gun_nation),
#   makeEstRDT2(out, treat=dates[6],name = 'Las Vegas',org= 'Gun Nation', y=gun_nation),
#   
#   makeEstRDT2(out, treat=dates[4],name = 'Atlanta + Boulder',org= 'Gun Store', y=gun_store)
# ) 
# write_csv(dat4meta, file='/Users/tyler/Dropbox/Mass Shootings and Tweets, Trends, and Donations/data/documentary/documentaries_meta.csv')

empty <- tibble(Coef = NA, lower = NA, upper = NA, shooting ='Sandy Hook',outcome = 'Gunned Down')
dat <- bind_rows(
  dat,empty
)

dat <- dat %>%
  mutate(shooting = factor(shooting, levels=rev(
    allshootings$shooting
  )))
dat$statsig <- factor(ifelse(dat$lower > 0,1,0))
dat$statsig[is.na(dat$statsig)] <- 0
dat$outcome[dat$outcome == 'Gunned Down'] <- 'Gunned Down (2014)'
dat$outcome[dat$outcome == 'NRA Under Fire'] <- 'NRA Under Fire (2020)'
dat$outcome[dat$outcome == 'Gun Nation'] <- 'Gun Nation (2016)'
dat$outcome[dat$outcome == 'Gun Store'] <- 'Gun Store (2019)'

dat$outcome <- factor(dat$outcome, levels=c(
  'Gunned Down (2014)','Gun Nation (2016)', 'Gun Store (2019)','NRA Under Fire (2020)'
))

calcMeta <- function(data, name, outcome, shooting){
  m.gen <- meta::metagen(TE = Coef,
                   seTE = se,
                   studlab = shooting,
                   data = data,
                   sm = "SMD",
                   fixed = FALSE,
                   random = TRUE)
  results <- tibble(shooting = name, 
                    Coef = m.gen$TE.random,
                    lower = m.gen$lower.random,
                    upper = m.gen$upper.random,
                    outcome = outcome,
                    statsig=factor(1))
  return(results)
}

meta1 <- calcMeta(dat %>% mutate(se = (upper - Coef)/1.96), 
         'Meta Estimate (Top 20)','Meta Analysis Estimate','Meta Analysis')
meta2 <- calcMeta(dat %>% mutate(se = (upper - Coef)/1.96) %>%
                   filter(shooting %in% top10), 
                 'Meta Estimate (Top 10)','Meta Analysis Estimate','Meta Analysis')
meta3 <- read_csv('datasets/lower_20_documentary.csv')
meta3$statsig <- factor(meta3$statsig)
meta3$shooting <- 'Meta Estimate (Bottom 20)' 
meta3$outcome <- 'Meta Analysis Estimate'

# now calculate meta analyses for race analysis
# race_est <- read_csv('datasets/shooting_race.csv')
# terciles <- quantile(race_est$mean.white, na.rm=T, c(0,0.33,0.66,1))
# race_est$tercile <- ifelse(race_est$mean.white < terciles[2],1,0)
# race_est$tercile <- ifelse(race_est$mean.white > terciles[2] & race_est$mean.white < terciles[3],2,race_est$tercile)
# race_est$tercile <- ifelse(race_est$mean.white > terciles[3],3,race_est$tercile)
# datnew <- left_join(dat, race_est, by='shooting')
# meta_top3 <- calcMeta(datnew %>% mutate(se = (upper - Coef)/1.96) %>% filter(tercile == 3),  
#                   'Meta Estimate (Top)','Meta Analysis Estimate','Meta Analysis')
# meta_top2 <- calcMeta(datnew %>% mutate(se = (upper - Coef)/1.96) %>% filter(tercile == 2), 
#                       'Meta Estimate (Middle)','Meta Analysis Estimate','Meta Analysis')
# meta_top1 <- calcMeta(datnew %>% mutate(se = (upper - Coef)/1.96) %>% filter(tercile == 1), 
#                       'Meta Estimate (Lower)','Meta Analysis Estimate','Meta Analysis')
# meta_tercile_out <- bind_rows(meta_top3, meta_top2, meta_top1) %>%
#   mutate(category = 'Documentaries')
# write_csv(meta_tercile_out, 'datsets/metaanalysis/meta_race_documentary.csv' )

plot1 <- dat %>%
  bind_rows(data.frame(
    Coef = NA,
    lower = NA, upper = NA, 
    shooting = na.omit(allshootings$shooting[!allshootings$shooting %in% unique(dat$shooting) & allshootings$include==1]),
    outcome = NA, statsig=NA
  ), meta1, meta2, meta3) %>%
  mutate(shooting = factor(shooting, levels=c('Meta Estimate (Bottom 20)',
                                              'Meta Estimate (Top 20)',
                                              'Meta Estimate (Top 10)',rev(allshootings$shooting)))) %>%
  ggplot(aes(shooting, y=Coef, ymin=lower, ymax=upper, shape=outcome, color=statsig)) +
  geom_hline(yintercept=0, linetype=2, color='red') +
  geom_errorbar(width=0, position=position_dodge(width=0.75)) +
  geom_point(position=position_dodge(width=0.75), size=2, fill='white') +
  coord_flip() +
  #facet_wrap (~outcome, ncol=1) +
  labs(title='',x='',y='Effect on Documentary\nWeb Streams')  +
  scale_color_manual(values=c('grey50','Black'), guide='none') +
  scale_shape_manual(values=c(21, 22, 23, 24, 16), name='', breaks=c(
    c('Gunned Down (2014)','Gun Nation (2016)','Gun Store (2019)','NRA Under Fire (2020)','Meta Analysis Estimate')
  )) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) 
plot1
ggsave(plot1, width=7, height=5, 
         filename = 'figures/documentaries_rddtools.png')

# low attention shootings

rm(list=ls())

pacman::p_load(lubridate, tidyverse, haven, gtrendsR, rdrobust, rddtools,
               broom, ggpubr, gridExtra)
setwd('put_your_wd_here')

# shootings
allshootings <- read_csv('datasets/list_of_shootings.csv',col_select = 1:4)
allshootings$date  <- mdy(allshootings$date)

# Frontline ----
df <- read_csv('datasets/frontline2.csv')
df$date <- ymd(df$date)
mindate <- min(df$date) + months(2)
maxdate <- max(df$date) - months(2)
dates <- allshootings %>%
  filter(date >= mindate & date <= maxdate & include==0)

# merge in guardian documentary
gua <- read_csv('datasets/gun_nation.csv') %>%
  dplyr::select(date, gun_nation)
df <- left_join(df, gua)

out <- readxl::read_excel('datasets/gun_store.xlsx', sheet=2)
out <- filter(out, Geography == 'US')
out$date <- ymd(out$Date)
out <- out %>% rename(gun_store = Views) %>%
  select(date, gun_store)
df <- full_join(df, out, by='date')

makeEstRDT <- function(out, treat, name, org, y, time=2){
  
  min <- treat - months(time)
  max <- treat + months(time) 
  out$running_day <- as.numeric(out$date - treat)
  trim <- out %>%
    filter(date < max & date > min)
  ggplot(trim, aes(running_day, gunned_web)) +
    geom_line() +
    geom_vline(xintercept=treat)
  y = trim %>% select(!!enquo(y)) %>% pull
  
  house_rdd <- rdd_data(y=y, x=trim$running_day, cutpoint=0)
  bw_ik <- rdd_bw_ik(house_rdd)
  reg_nonpara <- rdd_reg_np(rdd_object=house_rdd, bw=bw_ik)
  fullest <- data.frame(Coef=reg_nonpara$coefficients,
                        lower=(reg_nonpara$coefficients - 1.96 * reg_nonpara$coefMat[2]),
                        upper=(reg_nonpara$coefficients + 1.96 * reg_nonpara$coefMat[2]),
                        shooting=name,
                        outcome=org)
  return(fullest)
}

makeEstRDT2 <- function(out, treat, name, org, y, time=2){
  
  min <- treat - months(time)
  max <- treat + months(time) 
  out$running_day <- as.numeric(out$date - treat)
  trim <- out %>%
    filter(date < max & date > min)
  
  y = trim %>% select(!!enquo(y)) %>% pull
  
  sd_est <- sd(y, na.rm=T)
  house_rdd <- rdd_data(y=y, x=trim$running_day, cutpoint=0)
  bw_ik <- rdd_bw_ik(house_rdd)
  reg_nonpara <- rdd_reg_np(rdd_object=house_rdd, bw=bw_ik)
  fullest <- data.frame(Coef=reg_nonpara$coefficients/sd_est,
                        lower=(reg_nonpara$coefficients - 1.96 * reg_nonpara$coefMat[2])/sd_est,
                        upper=(reg_nonpara$coefficients + 1.96 * reg_nonpara$coefMat[2])/sd_est,
                        shooting=name,
                        outcome=org)
  return(fullest)
}

gunned_down <- vector('list', nrow(dates))
for(i in 1:nrow(dates)){
  gunned_down[[i]] <- makeEstRDT(df, treat=dates$date[i],name = dates$shooting[i], org= 'Gunned Down', y=gunned_web)
}

nra_web <- vector('list', nrow(dates))
for(i in 1:nrow(dates)){
  nra_web[[i]] <- tryCatch({makeEstRDT(df, treat=dates$date[i],name = dates$shooting[i], org= 'NRA Under Fire', y=nra_web)},
                           error=function(e){NA})
}

gun_nation <- vector('list', nrow(dates))
for(i in 1:nrow(dates)){
  gun_nation[[i]] <- tryCatch({makeEstRDT(df, treat=dates$date[i],name = dates$shooting[i], org= 'Gun Nation', y=gun_nation)},
                              error=function(e){NA})
}

gun_store <- vector('list', nrow(dates))
for(i in 1:nrow(dates)){
  gun_store[[i]] <- tryCatch({makeEstRDT(df, treat=dates$date[i],name = dates$shooting[i], org= 'Gun Store', y=gun_store)},
                             error=function(e){NA})
}

dat <- bind_rows(
  bind_rows(gunned_down%>% discard(is.na(.))),
  bind_rows(nra_web %>% discard(is.na(.))), 
  bind_rows(gun_nation%>% discard(is.na(.))), 
  bind_rows(gun_store%>% discard(is.na(.)))
) %>% na.omit() %>%
  remove_rownames()

dat <- dat %>%
  mutate(shooting = factor(shooting, levels=rev(
    allshootings$shooting
  )))
dat$statsig <- factor(ifelse(dat$lower > 0,1,0))
dat$statsig[is.na(dat$statsig)] <- 0
dat$outcome[dat$outcome == 'Gunned Down'] <- 'Gunned Down (2014)'
dat$outcome[dat$outcome == 'NRA Under Fire'] <- 'NRA Under Fire (2020)'
dat$outcome[dat$outcome == 'Gun Nation'] <- 'Gun Nation (2016)'
dat$outcome[dat$outcome == 'Gun Store'] <- 'Gun Store (2019)'

dat$outcome <- factor(dat$outcome, levels=c(
  'Gunned Down (2014)','Gun Nation (2016)', 'Gun Store (2019)','NRA Under Fire (2020)'
))

calcMeta <- function(data, name, outcome, shooting){
  m.gen <- meta::metagen(TE = Coef,
                   seTE = se,
                   studlab = shooting,
                   data = data,
                   sm = "SMD",
                   fixed = FALSE,
                   random = TRUE)
  results <- tibble(shooting = name, 
                    Coef = m.gen$TE.random,
                    lower = m.gen$lower.random,
                    upper = m.gen$upper.random,
                    outcome = outcome,
                    statsig=factor(1))
  return(results)
}

meta1 <- calcMeta(dat %>% mutate(se = (upper - Coef)/1.96), 
                  'Meta Analysis Estimate','Meta Estimate (All)','Meta Analysis')
meta_save <- calcMeta(dat %>% mutate(se = (upper - Coef)/1.96), 
                  'Meta Analysis Estimate','Meta Estimate (Bottom 20)','Meta Analysis')

write_csv(meta_save, file = 'datasets/lower_20_documentary.csv')

plot1 <- dat %>%
  bind_rows(data.frame(
    Coef = NA,
    lower = NA, upper = NA, 
    shooting = na.omit(allshootings$shooting[!allshootings$shooting %in% unique(dat$shooting) & allshootings$include==0]),
    outcome = NA, statsig=NA
  ), meta1) %>%
  mutate(shooting = factor(shooting, levels=c('Meta Analysis Estimate',rev(allshootings$shooting)))) %>%
  ggplot(aes(shooting, y=Coef, ymin=lower, ymax=upper, shape=outcome, color=statsig)) +
  geom_hline(yintercept=0, linetype=2, color='red') +
  geom_errorbar(width=0, position=position_dodge(width=0.75)) +
  geom_point(position=position_dodge(width=0.75), size=2, fill='white') +
  coord_flip() +
  labs(title='(F) Documentary Streams',x='',y='Effect on Documentary\nWeb Streams')  +
  scale_color_manual(values=c('grey50','Black'), guide='none') +
  scale_shape_manual(values=c(21, 22, 23, 24, 17, 20), name='', breaks=c(
    c('Gunned Down (2014)','Gun Nation (2016)','Gun Store (2019)','NRA Under Fire (2020)','Meta Estimate (Top 10)','Meta Estimate (All)')
  )) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) 
plot1
plot1 %>%
  ggsave(width=7, height=5, filename = 'figures/documentaries_rddtools_low_attn.png')
